Flow of power-law fluids in self-affine fracture channels 
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The two-dimensional pressure driven flow of non-Newtonian power-law fluids in self-afline fracture 
channels at finite Reynolds number is calculated. The channels have constant mean aperture and 
two values ("=0.5 and 0.8 of the Hurst exponent are considered. The calculation is based on the 
lattice-Boltzmann method, using a novel method to obtain a power-law variation in viscosity, and 
the behavior of shear-thinning, Newtonian and shear-thickening liquids is compared. Local aspects 
of the flow fields, such as maximum velocity and pressure fluctuations, were studied, and the non- 
Newtonian fluids were compared to the (previously-studied) newtonian case. The permeability 
results may be collapsed into a master curve of friction factor vs. Reynolds number using a scaling 
similar to that employed for porous media flow, and exhibits a transition from a linear regime to a 
more rapid variation at Re increases. 



I. INTRODUCTION 

An understanding of flow and transport processes in 
geologically disordered media is necessary for the effi- 
cient extraction of fluids from underground hydrocarbon 
reservoirs. Situations where flow proceeds through net- 
works of connected fractures are particularly attractive, 
because the throughput is generally much higher than 
may be achieved through intergranular porosity alone 
[H> S H) 13 • An important feature of subsurface frac- 
tures, which considerably complicates the problem, is 
that the surfaces of naturally fractures rocks are not 
smooth or even randomly rough, but rather are highly 
correlated self-affine fractals [5|. A second complication 
in the analysis is that typical reservoir fluids are often 
complicated mixtures, which exhibit non-Newtonian flow 
behaviors such as shear-thinning or shear-thickening. Yet 
a third difficulty is that the subsurface fracture flow often 
involves much higher velocities than in the intergranu- 
lar case, and the common simplification of low-Reynolds 
number linear flow is inapplicable. 

In this paper we use lattice Boltzmann calculations to 
elucidate the combined effects of self- affinity, non-linear 
rheology and finite inertia in fluid flow through a sin- 
gle fracture. Previous authors have considered subsets 
of these complications, but not all three simultaneously. 
The flow of Newtonian fluids in self-affine fractures at 
both low [fj, |7[ and finite [8J Re has an extensive litera- 
ture. Some controlled experiments on shear-thinning flu- 
ids in self-affine fractures at low Re have been reported 
9]. Lastly, experiments and phenomenological models 
for non-linear fluid motion in intergranular porous me- 
dia at various Re are available We anticipate that 
flow in a fracture can be characterized in a manner sim- 
ilar to the latter problem, since in both cases the key 
effect is that the random solid boundary of the flow do- 
main causes streamlines to wind around. One simplifica- 
tion which we can exploit, however, is to focus on two- 
dimensional flows. It is well known that the flow of a 
single fluid in a straight channel differs only in detail be- 
tween two and three-dimensional cases, and furthermore, 



in porous media flow in the analogous intergranular case, 
one sees the same flow laws for both two and three di- 
mensional geometries. 

The approach taken in the paper follows the lines of 
our previous studies of permeability [fj and transport 
in self-affine fractures based on the lattice-Boltzmann 
method, along with a procedure for incorporating power- 
law viscosity variation similar to that developed previ- 
ously [ll[ . The discussion of inertial effects is influenced 
by previous studies for the case of a Newtonian fluid in 
intergranular porous media @. The fracture surface is 
generated numerically by a Fourier transform algorithm 
and discretized on the regular lattice used in the flow 
problem. The upper and lower fracture surfaces bound 
the allowed nodes in the flow domain, a bounce-back con- 
dition enforces the no-slip boundary condition, and con- 
stant forcing provides a pressure driven flow. For power- 
law fluids, the lattice-Boltzmann relaxation time is ad- 
justed locally in space and time to provide the desired 
relation between stress and strain. The relation between 
imposed pressure drop and total fluid flux provides the 
permeability, and the local flow fields are analyzed to 
discuss the velocity, pressure and shear stress variations. 
Some background on the flow geometry and calculational 
method is presented in Section|TTl the local analysis of the 
flow fields in Sec. IIII1 the discussion of permeability is in 
Sec. IIV1 and we summarize in Sec. W\ 



II. BACKGROUND 

A. Self-affine roughness 

In this subsection we review the characterization of 
self-affine fractures and their numerical implementation. 
We consider a fracture surface without overhangs, i.e., 
the surface height h(x, y) is a single- valued function of the 
two coordinates r = (x,y) lying in the mean plane of the 
surface. A self-affine fractal surface is one which displays 
different scaling along the different spatial directions [12], 
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a statistical self-similarity under the transformation 

x — ► \x and y — > Ay => ft-(r) — > A^/i(r) (1) 

where £ is the Hurst or roughness exponent. Observa- 
tions of a variety of naturally fractured rock surfaces in 
different fracture modes yield just two common values of 
£, approximately 0.5 and 0.8. We further assume that 
the surface has spatial isotropy in its mean plane. The 
surface is further characterized by the amplitude of the 
roughness, or equivalently the prefactor Co in the height- 
height correlation function, 

<[Mr + A)-Mr)] 2 HCo(|A|/^) 2 < (2) 

where the intrinsic length scale £ might be the grain 
size in experiment or the lattice spacing in a calculation. 
In practice we generate self-affine surfaces using the a 
Fourier synthesis method [l3[ as in 0. 

A self-affine fracture channel is made of two comple- 
mentary self-affine surfaces separated by a gap, and in 
some cases the surfaces are shifted relative to each other 
parallel to the mean plane. The statistical properties 
of the fracture are specified by the Hurst exponent, the 
mean aperture between two surfaces, the shift distance, if 
any, and by the amplitude of the roughness. The height 
fluctuations of a single self-affine surface increase with 
its lateral extent L, so that the difference between the 
maximum and minimum heights scales as {L/ip, and 
we consider the limit H <C R < L, as shown for a typical 
fracture in Fig. [TJ 




FIG. 1: Geometry of a typical self-affine fracture composed 
of two complimentary self-affine surfaces with £=0.8. 

Note that the effective flow diameter of the fracture 
varies along its length and can be much smaller than 
the mean aperture, due to the tortuosity of the channel. 
When a lateral shift is present, the aperture varies locally 
as well, and furthermore if H is too small the sides of the 
fracture may overlap. 

B. The lattice-Boltzmann method 

Since the flow domain is bounded by highly irregular 
surfaces, the lattice Boltzmann method |14| is particu- 
larly convenient for fluid mechanical calculations, since 
the excluded solid region may be simply specified by 
a mask. If fi (x, t) is the velocity distribution function 
(VDF) for particles moving in direction i at lattice site 



x at time t, then the discrete Boltzmann equation which 
evolves the distribution is 

/,(x + e h t + 1) = /i(x, t) + ^(/(x, *)), (3) 

Here the e, are unit lattice vectors, the lattice spacing 
and the time step are both set equal to one, f^(/(x, <)) 
is collision operator which redistributes the VDF along 
different directions, and the spatial and temporal step 
discretizes in single unit. To recover the Navier-Stokes 
equation of fluid flow starting from Boltzmann equation, 
moments of the VDF satisfy the constraints 

p = Eifi pu = ^2f l e l <x = -pc£l-(l-— )^/ie»ei 

i i 

which relate the distribution function to the continuum 
density, velocity and stress fields, and where c s is the 
sound speed. The collision operator is treated in the 
BGK approximation using a single characteristic relax- 
ation time r, 

^(/(x,t)) = -i(/ i (x,t)-/f(x,t)), (5) 

T 

where /^ 9 (x, t) is equilibrium distribution function. The 
relaxation time r is related to the kinematic viscosity of 
the fluid by v — (2r — l)/6. To simulate a constant pres- 
sure gradient we add a a constant body force term to the 
right hand side of eqO equation eqJ3] and [5] to Navier- 
Stokes equation. More details may be found in [l4j], and 
a recent review of flow simulations in this context is pre- 
sented by Verberg and Laddjlq] . 



C. Power- law fluids 

The basic idea in extending the lattice Boltzmann 
method to power-law fluids was presented by Aharonov 
and Rothman and consists of adjusting the relax- 
ation time r locally so as to achieve the desired ratio of 
stress to strain rate. Here we consider power-law fluids 
using a generalized Newtonian model, as in (llj . where 
the relation between the stress tensor a a p and the strain 
rate tensor D a p = l/2((d0U a +d a up) is similar to that for 
Newtonian fluids, a a p = 2(j,D a p, but the local viscosity 
p is a function of the invariants of the strain rate tensor. 
We consider power-law fluids, /j = m A f n , where the 
case < n < 1 corresponds to shear-thinning, n > I cor- 
responds to shear-thickening, and n — 1 recovers linear 

Newtonian fluids, where the local shear rate 7 is related 

1/2 

to the second invariant of Dij via 7 = (2D : D) 1/z . The 
procedure in [TTI ] was to obtain the strain rate tensor by 
numerical differentiation of the previously calculated ve- 
locity field, then determine the appropriate local viscos- 
ity and thence the local relaxation time. Here we adopt 
a different procedure: in the lattice Boltzmann method 
the strain rate tensor is directly related to the velocity 
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distribution function by [l7j 



which should in turn equal <7 Qj g/2/i, there is a constraint 
on the fi which is solved by iteration. 

To validate the formulation of power-law fluids given 
above, we calculate the velocity profile for pressure- 
driven flow in a smooth-walled channel of constant aper- 
ture (a Hele-Shaw cell), which may be compared to an 
analytic solution of the Navier-Stokes equation. Apply- 
ing a pressure gradient AP/L = —G in the x-direction, 
the velocity for a power-law fluid with rheological param- 
eters 77i, n > as above in a channel of width H is 
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(7) 

We also record the mean velocity u and the fluid flux 
Q (per unit length in the passive third direction), which 
will be useful below: 

(8) 

In the simulation, we begin with zero velocity and in- 
tegrate Eq. (J3J to steady state, using the convergence 
criterion 



E 



||u(x, t) — ti(x, t — 1) 
||«(x,t)|| 



< 1.0 x 10" 



(9) 



For power-law indices n=0.75, 1.0 and 1.25, rn=0.01, 
and pressure gradient G = 1 x 10 -6 we obtain the profiles 
shown in Fig. [2l which agree with theory. In practice, as 
with any numerical method, computational instabilities 
may occur for substantially different values of the pres- 
sure gradient and fluid index, but the algorithm could be 
extended there using techniques such as multi-time step 
relaxation for the local shear viscosity [l8[ . 



III. LOCAL ANALYSIS OF THE FLOW FIELD 

We wish to examine how the local flow behavior varies 
with the rheology of the fluid, at different geometrical 
features of a self-affine channel. We focus on a single 
realization of the fracture, shown in Fig. [TJ and vary the 
power law index n and the pressure gradient G. The 
complete simulation box has length L — 256 in the flow 
direction and width W = 80, in terms of the (unit) lattice 
spacing, and the (constant) vertical aperture is H = 20. 
A uniform pressure gradient is applied everywhere along 
the channel, as above, and periodic boundadry conditions 
are applied in the flow direction. Local minimum in the 
effective width (normal to the average flow) occur around 
a;=55, 110, and 240 where mass conservation implies the 
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FIG. 2: Velocity profiles of power-law fluids with m = 
0.01, n = 0.75,1.0,1.25 in a Hele-Shaw cell with pressure 
gradient G = 1 x 10~ 6 . The points are simulation re- 
sults while the solid lines are the analytical solution in (01 • 
The maximum velocities for the three fluids 
0.006, 0.048,0.169, respectively. 



velocity magnitude will be a maximum, irrespective of 
the rheology of the fluid. In Fig. [31 we show velocity fields 
and streamlines for the three fluids along a segment of the 
fracture channel 20 < x < 100 in Fig. [1] which includes a 
constriction, for applied pressure gradient G = 1 x 10~ 6 . 
As we see, the streamlines are tortuous and very roughly 
follow the channel walls, although recirculating eddies 
(closed vortices) may occur where the channel exhibits 
side branches or dead-end regions. Indeed, at the present 
flow rate an eddy appears in the shear-thickening case 
but not the others, presumably because the velocites are 
higher in that case. 



A. Velocity field 

First we examine the variation of maximum absolute 
velocity along the channel, in order to show how the fluid 
rheology influences the earlier results of Skjetne et al. 
Q for the Newtonian case. More precisely, for each x 
along the channel we compute the maximum over y of 
|u(x, y)\, although we would have reached the same qual- 
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FIG. 3: Segment of Velocity vector field with stream- 
lines of the flow for power-law fluid with m = 0.01, n = 
0.75(top), 1.0(middle), 1.25(bottom) and the pressure gradient 
applied is G = 1 x 10~ 6 . The segment extends from x = 20 
to x = 100. 



itative conclusions had we considered the maximum over 
y of u x (x, y). Calculations were performed for three val- 
ues of the pressure gradient, G = 1 x 10~ 6 , 5 x 1CP 5 
and 2 x 10~ 4 , which correspond to Reynolds numbers 
Re — 0.95, 37.0 and 92.7, respectively, for the Newto- 
nian fluid. Since the viscosity varies within the channel 
for the shear thinning and thickening fluids, there is no 
unique definition of Re in those cases, although a conve- 
nient choice will be introduced in Section IIVI for scaling 
purposes. 

The resulting plots of maximum velocity are shown in 
Fig. [H where each velocity is normalized by the aver- 
age streamwise flow velocity (referred to as the inter- 
stitial velocity u* in ■$]), which equals the flux divided 
by the channel width. Obvious peaks appear at the po- 
sitions of the visible constrictions in the channel near 
x = 55, 110 and 240, reflecting the narrowed aperture 
there. The normalized peak heights are fairly insensitive 
to the Reynolds number, although away from the peaks 
the trend is for maximum velocity to increases with Re. 
Note that for a flat channel, the the normalized maxi- 
mum absolute velocity would equal 1.5, so the values of 
5 or more seen here are a substantial enhancement. The 



peaks are not all closely correlated with channel constric- 
tions, however: near x = 70 and 130 maximum velocity 
peaks occur, but at these locations the channel is expand- 
ing just downstram of a constriction. It is also possible 
to calculate a "maximum velocity trajectory" , following 
Q, as the set of (x, y) gridpoints which at each x has the 
y-value corresponding to the position where the maxi- 
mum velocity occurs. For the most part our observations 
concerning the behavior of these trajectories is similar to 
that reported in this reference, but we do not observe 
the line-length of this trajectory decreasing monotoni- 
cally with Re. 
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FIG. 4: Maximum absolute velocity along the fracture chan- 
nel for shear- thinning(top, n=0.75),Newtonian(mi<i(i/e, n = 
1.0), and shear- thickening( fcottom, n=1.25) fluids for various 
applied pressure gradient G. Each maximum velocity curve is 
normalized by the corresponding the average flow velocity 
in the a;-direction. 

Comparing the other fluids to the Newtonian case, we 
see in Fig. |4] that the global maximum absolute velocity 
always occurs at the narrowest constriction near x = 110 
and the other primary peaks always occur at the same 
positions, x = 55 and 240, as well. Furthermore, each 
peak has roughly the same (normalized) velocity value. 
In the shear-thinning case, both the variation in x away 
from the peaks/constrictions and the variation with pres- 
sure gradient are weaker than in the other cases, which 
may be attributed to the fact that typical velocities in 
the fracture are smaller in this case, and inertial effects 
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play a weaker role. In the shear-thickening case, where 
typical velocites are larger, the maximum velocity values 
are larger off the peaks values, and furthermore exhibits 
rather more variation with x and Re than the other flu- 
ids. 

The probablity distribution of velocity magnitudes is 
also of interest @ , since the presence of low and high ve- 
locity components strongly influences mixing processes 
and transport of passive tracers and suspended particles 
[l9| . Histograms of the observed absolute value of the 
velocity for the three fluids at various pressure gradients 
are shown in Fig. [5j In all cases there is a peak near the 
origin, which reflects the numerous low-velocity zones in 
the crevasses at the fractures walls, along with a higher- 
velocity peak resulting from the rapid flow in the channel 
constrictions. The latter moves out to higher values as 
the pressure gradient increases (note the normalization 
by in the figure) Once again, the shear-thickening 
case behaves somewhat differently than the other two flu- 
ids, showing a less prominent and broader "constriction 
peak," and more variation with G. 



B. Pressure and stress field 

The distribution of pressure and stress in the fluid 
are important for non-Newtonian rheology, and in con- 
sidering possible erosive processes on the fracture walls. 
To contrast the behavior of the different fluids, Fig. [5] 
shows the pressure "fluctuations" along the channel for 
the three power-law fluids n = 0.75, 1.0 and 1.25. The 
fluctuation p' is the deviation in pressure from the im- 
posed linear gradient, which would vanish identically in 
a Hele Shaw geometry. In the figure, the fluctuation 
has been normalized by the imposed pressure difference, 
Ap = GL, and averaged over the channel width. For 
all three fluids, the pressure fluctuation are most signifi- 
cant in the vicinity of the main constrictions in the chan- 
nel where the fluid accelerates, rising just before each 
constriction's location and dropping rapidly as it is tra- 
versed. Some additional structure arises at positions x of 
bends in the flow path, another source of fluid accelera- 
tion. Again, the shear-thinning and Newtonian fluids be- 
have somewhat similarly, while the variation is strongest 
in the shear-thickening case. 
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FIG. 5: Distribution of normalized absolute velocity in the 
whole self-afnne fracture flow domain for different power-law 
fluids with pressure gradients G = l.Oe — 6, 5.0e — 6 and 
l.Oe — 5 (top to bottom). 
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FIG. 6: Normalized pressure fluctuation along the channel for 
different power-law fluids, at pressure gradient G = l.Oe — 5. 



The variation in fluctuation with imposed gradient is 
shown in Fig. and indicates the expected general in- 
crease in magnitude with G along the channel. 

To assess the effects of the flow on the fracture wall, 
we first calculate the average force exerted by the fluid 
on the wall, 



F = 



d£n • <t 



(10) 



where the integral runs over the fracture surface (a curve 
in this two-dimensional calculation), and n is the local 
normal to the wall. The force is then decomposed into 
x and y components, representing the average drag and 
thrust on the wall, respectively, and then normalized by a 
typical inertial pressure pu 2 /2 times the nominal surface 
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FIG. 7: Pressure fluctuations of different fluids with power 
n = 0.75, 1.0, 1.25 and m = 0.01 along a self-affine fracture 
channel under different applied pressure gradients. 



area L X 1, to give drag and thrust coefficients 



F, 



Lpv?/2 



Lpu z /2 



(11) 



Note that aside from the (reasonable) use of the iner- 
tial pressure, the remainder of the normalization is some- 
what arbitrary but a fixed constant for each fracture, and 
mainly serves to provide dimensionless drag and thrust 
coefficients. The drag and thrust forces for the lower and 
upper walls of the channel are similar but not identical 
because of the asymmetry of the fluid-solid boundary, 
and for definiteness we present only the forces on the 
lower wall. 

The results of calculating the drag and thrust coeffi- 
cients is shown in Figs. [5] and [5] for the three fluids with 
exponents n = 0.75, 1.0, 1.25. In all cases, the coefficients 
exhibit simple power-law behavior, provided G is not too 
large, and the transition to a different behavior at larger 
G may be associated with the onset of inertial effects (see 
the following section). This form of scaling behavior re- 
sult is consistent with the experimental results reported 
in [l(| , and the values of the slopes found in the log- log 
plots in the low-G range, -1.67, -1.02, -0.62 for n = 0.75, 
1.0 and 1.25, respectively, for both drag and thrust, may 



be understood from the following argument. 

If inertial effects are absent, one expects the scaling be- 
havior in a rough channel to be the same as in a straight 
channel. In that case, from Eq. ([7]) one has u ~ G 1 /™, 
and therefore Vw ~ G 1 /" as well, so that /i ~ |Vu| n_1 ~ 
Q(n-i)/n_ rpjjg (j ra g an d thrust forces are proportional 
to the stress, a ~ /zVu ~ G ( " _1)/n+1/n - G 1 . The drag 
and thrust coefficients are then d,t ~ F x _ y /u 2 ~ a/u 2 ~ 
giving exponents -5/3, -1 and -3/5, respectively, 
for the three fluids. 
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FIG. 8: Drag factor d for power-law fluids in a self-affine 
fracture channel as a function of applied pressure gradient. 
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FIG. 9: Thrust factor t for power-law fluids in a self-affine 
fracture channel as a function of applied pressure gradient. 



IV. PERMEABILITY 

Next we consider global behavior - the permeability of 
a self-affine fracture channel. Our discussion is colored 
by analogies to flow in intergranular porous media, so we 
first recall the situation in that system [20j . For Newto- 
nian fluids in the low Reynolds number limit, the defi- 
nition of intergranular permeability is given by Darcy's 
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law, (u) = —(k/fi)Vp, where (u) is the average flow ve- 
locity and p the average pressure. The average in ques- 
tion could be a volume average or an ensemble average, 
and for a flow which is macrosopically unidirectional, an 
operational definition of permeability is k = fiQL/AAp 
where Q is the flux through a sample of cross-sectional 
area A and length L. In a two-dimensional situation, the 
area is replaced by the width W, and Q is the flow per 
unit length in the third direction. A definition identi- 
cal to the latter case may be used for the permeability 
of low Reynolds number Newtonian flow in a fracture. 
Both finite Reynolds number flow and non-Newtonian 
fluid rheology modify this description. We first consider 
the effects of inertia, and then examine how permability 
relates to the fracture morphology. 



A. Inertial effects 

At higher flow rates when inertial effects appear, the 
relation between pressure difference and average velocity 
or flux becomes nonlinear and one may write 



Ap = aQ+[/3Q 2 or 7 Q 3 



(12) 



where a incorporates the Darcy permeability, and the 
term in brackets is the inertial correction, with a, f3 > 0. 
At high Q the quadratic or "Forchheimer" term applies, 
but in the transitional region where the Reynolds number 
is small but finite, a cubic dependence is found. This pic- 
ture is supported by experiments, analytic calculations, 
and numerical simulations plj . 

The flow of a Newtonian fluid in a self-affine fracture 
can be described in identical terms, as shown by the nu- 
merical simulations of Skjetne et al. 8] which exhibit 
the same transitions between flow regimes indicated in 
Eq. (fl2|) . In extending the discussion to power-law flu- 
ids, the first issue is to choose the appropriate power 
of Q. The exact solutions for Hele Shaw flow given in 
Eq. ([8|) have the scaling behavior G ~ Q n where G is 
the applied pressure gradient (the relevant pressure for 
macroscopic behavior) and n the power-law index. In a 
rough fracture, one would naturally expect an identical 
relation, albeit with a modified coefficient, at low G, and 
then at larger G inertial effects would be expected to pro- 
duce (positive) terms involving higher powers of Q. To 
test this idea, note that we are concerned here with the 
statistical behavior of self-affine fractures, rather than 
the details of flow in one particular geometry which was 
relevant in the previous section, so an ensemble average 
over six realizations of the fracture surface is used. The 
simulation results are shown in Fig. [TO] and indeed show 
a G ~ Q n scaling behavior at low G. The Newtonian 
n=l plot shows this behavior clearly since G/Q n is con- 
stant at low Q, whereas in the other cases, the expected 
behavior is present at sufficently small Q as indicated 
in the alternative plots in the insets of G vs. Q. The 
need for different plotting variables arises becuase in the 
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FIG. 10: Relation between imposed pressure gradient and 
fluid flux for power law fluids: n = 0.75, 1.0 and 1.25 (top to 
bottom). 



non-Newtonian cases, the flow rate fluctuates substan- 
tially at low G and division by Q n is numerically un- 
stable. Beyond the quasi-linear regime, the Newtonian 
case shows the expected trasition to a Forchheimer flow 
regime G ~ Q 2 at larger forcing, and the shear-tickening 
fluid shows a somewhat analogous behavior G ~ Q 2n . 
The shear-thinning fluid is not described by a simple 
power law at large G, and we are not aware of any the- 
oretical treatment of this problem, so we simply report 
the numerical results. 

To understand the numerical coefficient in the flow re- 
sults, the fracture-modified Darcy permeability, we again 
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fluid index fco k 
n = 0.75 2.99 0.373 
n = 1.0 33.2 4.73 
n = 1.25 142 20.6 



TABLE I: Effect of roughness and tortuosity on the low 
Reynolds number permeabilty: fco and fc are the permeabili- 
tites (defined in Eq. (|13p for a Hele Shaw cell and a self-affine 
fracture of the same mean aperture, respectively. 



how such ad hoc factors might be interjected here, so in- 
stead we collapse the data using a simple constant factor 
which varies from fluid to fluid, and the result is shown 
in Fig. [TTJ Two different values of the Hurst exponent 
are shown, and in both cases we see anF~ 1 /Re scaling 
at low Re, a transition at Re ~ 1 — 10 and perhaps a 
constant friction factor at larger Re. Unfortunately, the 
calculations cannot be extended into the latter regime 
using the present method (a particluar implementation 
of the lattice Bolzmann technique) becuase numerical in- 
stabilites arise. 



refer to the Hele Shaw case and define 



A; = 



m 1 ' n u 
G 1 /™ ' 



(13) 



Since the roughness and tortuosity of the fracture cause 
the stramlines to bend and viscous dissipation to in- 
crease, the permeability should be reduced compared to 
a smooth and flat Hele Shaw geometry of the same aper- 
ture. In Table I, the various permeabilites are compared, 
and a reduction by a factor 6-7 is found. 

So far, we have expressed the pressure gradient G in 
terms of the flux Q, because these quantities are well de- 
fined in the present simulations. However, for general 
purposes It is preferable to use a dimcnsionless quan- 
tity such as the Reynolds number as the independent 
variable, but the definition of Re for power-law fluids is 
not entirely obvious for power-law fluids because the the 
viscosity varies over the flow domain. One way to com- 
bine the results for different fluids is based on an anal- 
ogy to the friction factor scaling laws for flow in pipes 
originally due to Nikaradze [22j], which can be extended 
to non-Newtonian fluids as shown by Metzner [23|. Re- 
call that for unidirectional flow of a Newtonian fluid of 
viscosity p in a pipe of diameter D, the mean velocity 
is u = GD 2 /32p and the shear stress at the wall is 
t w = GD/A, so if one defines the conventional friction 
factor as / = T w /±pu 2 , then one finds / = 16/i?e where 
Re — puD / fx. Experiments follow this scaling law up to 
a value of Re that depends on the roughness of the pipe, 
and at larger values of Re, f levels off. An analogous cal- 
culation for Hele Shaw flow using the aperture H instead 
of the diameter D gives r w — GH/2 and / = 12/ Re. The 
power-law generalization is to use the latter form for t w , 
along with Eq. f8j to express the pressure gradient in 
terms of the mean velocity, and yields 



12 
Re 



if Re 



6pu 2 - n H n /m', 



(14) 



where m' = m(2(2n + l)/n) n . This choice of variables 
is not the last word, because in the analogous intergan- 
ular porous medium case where a similar approach has 
been taken [Io( ]. extra constant factors such as functions 
of the porosity or the "dynamic specific surface area" are 
introduced into the friction factor and Reynolds num- 
ber definitions to promote data collapse. It is not clear 
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FIG. 11: Friction factor of self-affine fracture channels of 
Hurst exponent ( = 0.5 and 0.8 as a function of Reynolds 
number defined as in [14] for power-law fluids with n — 
0.75,1.0,1.25. 



B. Morphology effects 

We now consider how the geometry of the fracture ef- 
fects the (low- Reynolds number) permeability for the var- 
ious fluids considered. First we investigate the effect of 
the Hurst exponent £ on the permeability, and to simplify 
the analysis we consider a fracture channel with one self- 
affine wall and one flat wall, as in 0] . For a fixed presure 
gradient G, we compute the flux as a function of the 
channel length L for the three fluids, and in Fig. [12l we 
first show the flow rate depletion (Qq — Q) /Qq vs. L for a 
fracture with ( = 0.8. Here Qq is the flux through a flat- 
walled channel fo the same mean aperture. Increasing the 
length allows for more fluctuation in the channel width 
(see Eq. (0) which increases the tortuosity and tends to 
decrease the flux. If the Hurst exponent of the channel's 
rough wall is instead ( = 0.5, the three fluids again be- 
have quite similarly, so it suffices to compare the behavior 
of different Hurst exponents for a single case, and in the 
lower panel of the figure we plot the flux depletions for 
the two exponents for the shear-thinning case. The fact 
that the flux depletion is greater for the £ = 0.8 chan- 
nel may be explained by noting that this exponent value 
corresponds to more fluctuation as a function of L than 
the 0.5 case, and therefore to a more tortuous channel. 
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sections. Using these relations in Eq[15]we have 
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FIG. 12: Flux variation with length in a channel with one 
self-affine and one flat wall for different fluids and Hurst expo- 
nents. The maximum aperture of the channel is J/ m ax = 64, 
and the applied pressure gradient is G — 1.0 x 10 -6 . Top: 
flux depletion for different fluids confined in a channel with 
£ = 0.8. Right: flux depletion for a shear-thinning fluid in 
channels of different f. 



To relate the flux to the channel aperture, we imagine 
dividing the channel into a sequence of nearly-straight 
sections, each of length £i, and writing the total pres- 
sure difference as the sum of the pressure drops in each 
section, using Eq. [8] for each. This reasoning yields 



n 

(15) 

where the summation is over the sections, and hi is the 
effective aperture and k is the length along the local flow 
direction in section i, and we have noted that Q is a the 
same in all sections. If 9i is the angle between the orien- 
tation of channel section i and the mean flow direction, 
then hi = HcosOi and = lj/cos0i, where H is the 
aperture and is the projected length of section i in 
the mean flow direction, assumed to be the same for all 



AP = 2mQ I" 



^(cos 0i)~ 



(2n+2) 



(16) 

This result generalizes Eq. 26 in [6[ to power-law fluids, 
and if we proceed as in that reference to evaluate the 
average over angles 9i we obtain 



Q - Q ~ Jj(2<-2)/C+(2n+l)/r. 



(17) 



where again Qq is the flux in a flat channel of the same 
aperture H. 

To test the relation [T3 we calculate the flow for frac- 
ture channels of length L = 256 with varying apertures 
H = 8, 12, 16, 20, 24, for fluid with m = 0.01, n=0.75, 1.0 
and 1.25, all at a pressure gradient AP/L = l.Oe — 6. 
Figure Q2] shows the flux depeltion (Qo — Q) as a func- 
tion of aperture. The points are the numerical results 
and the solid lines are fitted curves, based on the ex- 
pected power-law exponents obtained from EafTTl which 
are 2.83, 2.5, and 2.3 for the three fluids. We see that 
the theoretical analysis is in excellent agreement with 
the data for the shear-thinning and Newtonian fluids 
(n=0.75 and 1.0), but the agreement is less statisfactory 
for the shear-thickening fluid, whose numerical exponent 
is closer to 2.5. A possible interpretation is that in the 
shear-thickening case, for the same pressure gradient the 
average velocity is larger than that for the other fluids, 
so that fluid inertia comes into play. 
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FIG. 13: Log-log plot of flow rate variations versus the 
aperture of self-affine channel for different fluids with power 
n = 0.75, 1.0 and 1.25. 



Finally, we consider an additional effect, a lateral shift 
between the two sides of a fracture, which might arise 
in practice due to geological processes. We begin with 
a fracture channel with complimentary sides and con- 
stant initial aperture H, and then shift one side along 
the mean plane by a distance d. The fracture aperture 
is now a function of position, H^ix), and effectively a 
spatial random function. We again compute the flux de- 
pletion relative to a flat channel having the same ini- 
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tial aperture, using six realizations of a self-affine frac- 
ture wall with Hurst exponent £ = 0.8. As shown in 
Fig. [14] the flux decreases somewhat faster than linearly 
with shift, by producing narrow gaps when proturbanccs 
on the two sides are brought closer to one another. The 
shear-thinning and Newtonian fluids have a fairly simi- 
lar behavior, while the reduction is twice as large in the 
shear-thickening case, perhaps again as a result of iner- 
tial effects. As in the previous discussion, using a differ- 
ent value £ = 0.5 for the Hurst exponent gives the same 
trends. 
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FIG. 14: Flow rate reduction due to lateral shift along the 
mean plane of fracture channel for different fluids with power 
n = 0.75, 1.0 and 1.25, and the Hurst exponent used here is 
C = 0.8. 



V. CONCLUSION 

Using a new implementation of the lattice Boltzmann 
method for power-law fluids, we have investigated their 
flow in two-dimensional self-affine fracture channels as a 
function of applied pressure gradient. Generally, fluids 
with different power-law index behave in a similar man- 
ner when their flow parameters are properly scaled, us- 
ing standard results for flow in constant-thickness chan- 
nels. Many previous results for Newtonian fluids in self- 
affince fractures are found to generalize in a straightfor- 
ward manner. However, shear-thickening fluids, which 
have higher velocities for the same pressure gradient than 
Newtonian or shear-thnning counterparts, are more sus- 
ceptible to inertial effects. 

With regard to the local flow fields, we first consid- 
ered the maximum absolute velocity as a function of dis- 



tance along the mean flow direction, which was found to 
fluctuates along the fracture channel due to its tortuos- 
ity and the variable effective aperture along the channel. 
The local maxima of this maximum absolute velocity oc- 
cur at points of narrowing or minimal effective aperture, 
and the range of maximum absolute velocity relative to 
the global mean velocity ranges from about 1.5 to 5.5. 
With increasing inertia, this normalized maximum abso- 
lute velocity increases for all power-law fluids to different 
degrees, with shear-thickening fluids having the largest 
effect and shear-thinning the least. As the pressure gra- 
dient increases, the normalized maximum velocities near 
the constrictions are relatively constant but outside these 
points velocities tend to increase. The variation in veloc- 
ity is greatest for a shear-thickening fluid and lest for 
shear-thinning. Pressure fluctuations along the channel 
increase with forcing for all fluids, and for a given pres- 
sure gradient increase with the power-law index n. 

The relationship between pressure gradient and flux is 
found to have the same functional form as for flow in a flat 
channel, Ap ~ Q n , when inertial effects are absent. At 
higher Ap, Newtonian fluids behave in the same way as 
in intergranular porous media, and shear-thinning fluids 
behave analogously, but the shear-thickening case does 
not show simple power-law behavior. It is possible to 
collapse all of the data on flux vs. pressure gradient into a 
universal friction factor curve. The variation of flux with 
system length was shown to scale with system length with 
an exponent algebraically related to the Hurst exponent, 
in a manner which generalizes the Newtonian case. 

The most interesting question raised by these results 
is the form of the flux-pressure gradient relationship in 
the regime of strong inertia in the non-Newtonian case. 
In this work, we were limited in the range of accessible 
pressure gradients by numerical instabiliites, and it is de- 
sirable to improve the algorithm so as to consider higher 
pressure gradients and further explore the dynamics of 
the inertial regime. An extension of these considerations 
to viscoelastic fluids is likewise highly desirable, but new 
ideas beyond the methods used here are needed. 
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